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We investigate the possible effects that deviations from kinetic equilibrium can have on massive 
particles as they decouple from chemical equilibrium. Different methods of solving the Boltzmann 
equation yield significantly different relic number densities of such particles. General considerations 
concerning the Dirac or Majorana structure of the particles are discussed. It is shown that non- 
equilibrium effects are small for particles decoupling while strongly non-relativistic, as will be the 
case for most cold dark matter candidates. 
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Deviations from thermodynamic equilibrium are of fundamental importance in the early universe. The examples 
are plentiful, for example the freeze-out of photon thermodynamic equilibrium due to hydrogen recombination Jl ] or 
the equivalent freeze-out of neutrinos slightly prior to electron-positron annihilation Q . These phenomena have been 
■ investigated many times in the literature, but almost always in the context of deviations from chemical equilibrium, 
that is, kinetic or scattering equilibrium is assumed to hold at all times. In the limit of Boltzmann statistics the single 
C*~) particle distribution then takes the form 

o ■ 

f{E) = e-< B -^/ T , (1) 

ON ■ where fi is a pseudo-chemical potential, independent of E. In recent years there has been some discussions as to the 
validity of this assumption, especially concerned with the decoupling of neutrinos in the early universe |3-M|. In this 
case it has been found by several authors that the effects are small, but not negligible. In the present paper we will 
i' discuss chemical versus scattering equilibrium in a general context, where we limit ourselves to generic expressions 
for the matrix elements. Dolgov has provided a general discussion of deviations from kinetic equilibrium caused 
by the overall expansion of the universe, looking only at elastic scattering. The case we treat is where both scattering 
and annihilation is taken into account to solve the full Boltzmann equation, however we resort to a numerical solution 
Ci • of the Boltzmann equation. 

In general, for Dirac fermions it can be assumed that s-wave (J = 0) annihilation is by far the most important 
for non-relativistic particles. The reason is the following: since particle and antiparticle are two different particles 
there are no anti-symmetry requirements and the annihilation can proceed via s-waves. Since higher order terms 
J > 1 are suppressed by factors of vf£ <C 1 [jl0|], where u re i is the relative velocity of the incoming particles, s-wave 
annihilation is dominant. For Majorana particles particle and antiparticle are identical objects and the combined 
wave function must be antisymmetric. Therefore annihilation proceeds primarily via p-waves (J = 1), meaning that 
the annihilation amplitude is suppressed relative to the Dirac case. In the case of s-wave annihilation the matrix 
element is not momentum dependent and the annihilation amplitude is constant. For p-wave annihilation there will, 
as mentioned, be a dependence on the square of the relative velocity of the two incoming particles. For non-relativistic 
particles we can therefore in general write the annihilation amplitude as 

£ |M i 2 = { j=i! ■ < 2 > 

where G is a dimensionless coupling constant. 

We shall use this expression for the matrix element even for relativistic decoupling, although that is of course at 
best only an approximation to the true matrix element. For example, in the case of massive Dirac neutrinos, the 
matrix element will be strongly energy dependent as long as the neutrinos are relativistic, but approach a constant 
as the neutrinos become non-relativistic JtJ. However, taking our approach, although quantitatively not accurate for 
relativistic decoupling, still gives a very clear picture of when non-equilibrium effects are important. 

The elastic scattering part will, in general, have a momentum dependent amplitude. This is for example the case 
for scattering of heavy neutrinos on massless particles via the weak interaction. Here, the squared matrix element is 
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proportional to the squared energy of the scatterer. For the sake of simplicity we shall deal only with three cases, the 
first being no elastic scattering at all and the second elastic scattering with a constant amplitude. The first case will 
represent an extreme upper limit to the size of non-equilibrium effects. The second case, albeit not realistic, serves 
to indicate how elastic scattering influences the shape of the massive particle spectra. The third case we shall treat 
is the case where full scattering equilibrium is maintained at all times. 



II. BOLTZMANN EQUATIONS 



When discussing non-equilibrium effects in the early universe the relevant equation is the Boltzmann collision 
equation. Since we are not discussing such issues as medium effects we shall use the single-particle Boltzmann 
equation. Also, we restrict the calculation to using Boltzmann statistics, which simplifies the calculations while 
leaving the essential physics intact. 

Then the Boltzmann equation takes the form [[yj 

df = C a „n + C ch (3) 

where 

H being the Hubble parameter H = R/R. On the right-hand side, C an n and C e \ represent annihilations and elastic 
scatterings respectively. These collision terms can be written generically as 

CUE/] = tt^t / d 3 p 2 d 3 p 3 d 3 p 4 (f 3 f 4 - hf 2 ) (5) 
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xS]T |Ml? 2 ^34<* 4 (Pi + P2 - P3 - Pi){2ir)\ 

as long as we are only concerned with 2-body collisions. Here we have d 3 p = d 3 p/ ((2ir) 3 2E). S is a symmetrisation 
factor of 1/2! for each pair of identical particles in initial or final states |1J], and ^2 \M\ 2 is the interaction matrix 
element squared and spin-summed, pi is the four-momentum of particle i. Now, the next issue is to actually perform 
the relevant collision integrals, ft is in general possible to integrate C an n and C e \ down to two dimensions, depending 
on the actual shape of ^ \M \ 2 . The details of these calculations are given in the Appendix. 

The usual approach to solving the above equation, Eq. (||), is to assume full scattering equilibrium (our third case). 
Then the elastic scattering term is identically zero at all times, and the Boltzmann can easily be integrated to yield 

H 

n + 3Hn — — (a\v\)(n 2 — nl q ), (6) 

where n is the actual number density and n cq is the equilibrium number density. is the velocity-averaged 

annihilation cross-section H|. 

As for the general expansion of the universe we shall assume that the universe is completely radiation dominated 
and further we will ignore entropy production in order to make the discussion completely generic. Including the 
correct expansion law by using the Friedmann equation and the equation of energy conservation will only perturb the 
general results slightly in most cases. Using the Friedmann equation |13| , 

H> - (7) 
one then arrives at the following time-temperature relation 



' MeV 

where g* = p/py, and p is the total radiation energy density. 



t(.s)-1.71. 9 ; 1/2 T Mo 2 v , (8) 



2 



III. NUMERICAL IMPLEMENTATION 



In order to investigate the non-equilibrium effects in a quantitative way we still have to calculate the remaining 
two dimensions of the collision integrals, Eq. (^), numerically. We have done this in order to estimate whether non- 
equilibrium effects can ever be important during particle freeze-out. Notice that we do not make any assumptions as 
to the size of the non-equilibrium effects. That is, we do not assume a priori that they are small, as was done for 
instance in Ref. 0] , but instead we perform a full numerical integration of the whole system. In all calculations we will 
assume that the massless particles are in full thermodynamic equilibrium because of their interactions with other light 
particles, an assumption which simplifies the calculations. For both Dirac and Majorana particles one can expect that 
the coupling constant will in general be the same for annihilation and scattering, e.g. G oc Gf for weak interactions 
fFor a non-relativistic Dirac neutrino, G = 32m 4 G F , assuming only the annihilation channel vh^h — * v e v e ) or G oc a 
for electromagnetic interactions. We can therefore scale the amplitude in units of a single coupling constant, G, so 
that 

J2\M\ 2 kG 2 . (9) 




FIG. 1. The final number density in comoving coordinates after complete decoupling. The upper panel shows the case where 
the interaction matrix element is constant, both for scattering and annihilation. The lower panel shows the case for a velocity 
dependent annihilation matrix element, still with a constant scattering matrix element. 

Clearly, the size of the non-equilibrium effects will depend on G. If G is very small, particles will decouple while still 
relativistic. This of course would not produce any non-equilibrium effects at all. On the other hand, if the particles 
are extremely non-relativistic before they decouple, the effects are also likely to be small for the following reason. For 
very non-relativistic particles in equilibrium, the number density is suppressed exponentially relative to the massless 
species. Therefore the annihilation rate is also suppressed exponentially, while the scattering rate is suppressed at 
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most by some power of p/m. This will increase the importance of elastic scattering relative to annihilation and lessen 
any non-equilibrium effects, even though the scattering on massless particles is essentially momentum conserving. A 
priori one therefore expects that non-equilibrium effects are the largest for semi-relativistic decoupling. 

In Fig. 1 we have plotted the final number density of massive particles in comoving coordinates as a function of a 
dimcnsionless interaction strength, defined as 



r = 5.20 x 10 20 m-lyg^^G 2 . 



(10) 



There is a pronounced difference here between the two cases of J — and J = 1. For strong interactions (large 
r), the final number density is much higher for the velocity dependent interaction because of the large suppression 
factor v 2 el /m 2 . On the other hand, for relativistic decoupling, the velocity dependent interaction is stronger by the 
same factor for the same value of T and therefore the final number density in this case is lower than for the constant 
matrix element. Also, the freeze-out process proceeds much faster in the J = 1 case than in the J = because the 
annihilation amplitude falls off with decreasing temperature. 




FIG. 2. Evolution of comoving number density y = nR 3 in a specific case, where we have taken F = 0.5. The full curve 
shows the constant matrix element (J = 0), the dotted shows the velocity dependent matrix element (J = 1), and the dashed 
line is the equilibrium number density. 



In Fig. 2 we show how the annihilation actually proceeds in the two cases. This example is taken for T = 0.5. 
Clearly, for the constant matrix-element decoupling is very prolonged, taking place over a wide range in T. For the 
velocity dependent matrix element, decoupling is much more sudden, essentially completing over a factor of two in T. 
This difference will of course give rise to very different non-equilibrium effects. 

Fig. 3 shows the deviation in final number density between relative to that found by using the integrated Boltzmann 
equation. For the solution of the full equation (the full lines) We see the expected trend, namely that relativistic 
decoupling leads to no significant non-equilibrium effects and the same being true for very non-relativistic decoupling. 
The effect has a maximum for semi-relativistic decoupling, where the difference between using the full Boltzmann 
equation and the integrated version can be as large as 10-15%. However, even for strongly non-relativistic decoupling 
there still is a non-negligible effect because of residual annihilations taking place after decoupling from both scattering 
and chemical equilibrium. Thus, one cannot entirely neglect non-equilibrium effects when doing this sort of calculation. 
Indeed, for the special case of massive neutrino decoupling it has been demonstrated previously in the literature that 
out-of-equilibrium effects can be significant ||-[|- Note also the interesting feature that for the case of a constant 
matrix element, the annihilation term will itself preserve scattering equilibrium if the massive particles are non- 
relativistic. This can be seen directly from the kernel in the annihilation integral, Eq. (A7). The dotted curve in the 
upper panel of Fig. 3 shows an example of this behaviour. Here, the size of the non-equilibrium effects also drops 



4 



for strong interactions because complete equilibrium is maintained until the particles are strongly non-relativistic, 
where after C an n preserves scattering equilibrium independent of the interaction strength. Fo r th e velocity dependent 
matrix element the situation is the reverse because the annihilation integral kernel, Eq. (|A9|), does not preserve 
kinetic equilibrium. Rather, high momentum states annihilate much faster than low momentum ones, so that kinetic 
equilibrium is quickly destroyed. This can lead to very large errors in the computed number density if one uses the 
approach of only taking into account the annihilation term in the collision integral. 
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FIG. 3. The deviation in number density depending on the calculational approach. The full line is the deviation found found 
by using the full Boltzmann equation relative to using the integrated equation. The dotted line is the difference between using 
the annihilation term only in the full Boltzmann equation and the integrated equation. 



The sign of the non-equilibrium effect is also opposite for the J = and J = 1 scenarios. For the constant matrix 
element low momentum states will, in general, annihilate faster than high momentum ones so that maintaining 
kinetic equilibrium leads to a lower annihilation rate. For the velocity dependent matrix element high momentum 
states deplete quickly so that annihilation is slowed unless they are refilled by scattering. This effect is known from 
the case of massive Majorana neutrinos || || . 

Finally, in Fig. 4 we show an example of the actual distribution functions obtained by different methods, taken 
as a snapshot at T/m = 0.15 for T = 0.5. What is shown is the ratio between the actual distribution function and 
the corresponding equilibrium distribution with the same number density. Obviously, the distribution one obtains 
by using the integrated equation is an equilibrium distribution per default. However there can be huge differences 
between that and the actual distribution function obtained from the full Boltzmann equation. Also, there is a large 
difference between the full solution and the one obtained by using only the annihilation term. 
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FIG. 4. Actual distribution function of the heavy particle. The full line corresponds to the solution of the full Boltzmann 
equation, the dotted to including only annihilation and the dashed to the equilibrium distribution in the integrated Boltzmann 
equation. The distributions have been taken at a temperature of T/m = 0.15 and for V — 0.5. 

IV. DISCUSSION 

Our results clearly show that for the case of semi-relativistic decoupling one cannot safely ignore the non-equilibrium 
effects pertaining to the deviations from scattering equilibrium. This will for instance be the case with MeV mass 
neutrinos which can influence Big Bang nucleosynthesis. 

However, for either relativistic or non-relativistic decoupling the effects are negligible. Thus, for most cold dark 
matter candidates which per definition decouple while strongly non-relativistic [], it is completely safe to assume full 
scattering equilibrium until long after decoupling from chemical equilibrium. 

As a simple example we shall consider a massive Dirac neutrino which is thermally produced in the early universe. 
This particle is an excellent CDM candidate if it has a mass in the 1 GeV region, around the Lee- Weinberg limit jl5| . 
The dimensionless interaction strength for this particle is given by 



2-3 g. 



-1/2 3 

m McV' 



(11) 



assuming that the only relevant annihilation process is vhvh v e v e . In Fig. 5 we have plotted F as a function of 
the heavy neutrino mass. For a mass of 1 GeV, T ~ 7 x 10 s <C 1, meaning that there are essentially no deviations 



An example of a CDM candidate with radically different properties is the axion which is produced in a condensate at T m, 
see for instance Ref. Ill 31 
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from kinetic equilibrium. As a note of caution we stress that we have assumed a scattering matrix element of the 
same magnitude as the annihilation matrix element which is incorrect. Rather, the scattering matrix element is 
2 I M | scattering — G 2 F m 2 T 2 which is lower by a factor of T 2 /m 2 . Even so there will be no discernible effect due 
to deviations from scattering equilibrium for such a heavy neutrino. The same will be true for super-symmetric 
dark matter particles which are likely to have a mass in the multi-GeV region or above [Tq| . However, there are 
other examples where non-equilibrium effects are sizeable, for example the previously mentioned example of an MeV 
neutrino. 
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FIG. 5. The dimensionless interaction strength F, defined in Eq. 
the mass. 



lOh, for a non-relativisitc Dirac neutrino as a function of 



It was also shown that taking into account only the annihilation term in the full Boltzmann equation can lead to 
large errors in the final number density if the matrix element is velocity dependent. This was for instance the case in 
the work by Kawasaki et al. ||. 

Finally, in the present paper we have not discussed possible deviations from equilibrium in the massless particle 
species induced by interactions with the massive particles. This effect is for instance of importance when the massless 
particles themselves are decoupling from equilibrium at the same temperature as the massive ones. An example is the 
decoupling of MeV r-neutrinos where the induced deviations from equilibrium in the electron neutrino distribution 
causes Big Bang Nucleosynthesis to change 0,0,^-^. 
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APPENDIX: COLLISION INTEGRALS 



This appendix deals with the phase-space integrations of the collision operator of Eq. (|5|) . The principal method 
of integration was developed in Ref. |F7| and in the following we only give a brief outline and further develop the 
integrations in the non-relativistic case. 

The collision integral in the general two-body collision integral can be reduced to two dimensions. Let the process 

be 
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Then the result is 



where 



1 + 2^3 + 4. (Al) 
CWACpi)] = J d P2 dp 3 [J(p 1 ,p 2 ,p 3 )}^ a 1 2- ] (f3h ~ A/a), (A2) 

± -27 -2pj-Q± 2 P2 (27 + pj + g| + g| + g) I 

a = , (A3) 

2pip 3 

1 = E 1 E 2 -E 1 E i -E 2 E^ (A4) 
Q = tn 2 + m 2 + — mf. (A5) 

The function J{pi,p 2 ,p 3 ) is the matrix element integrated over all variables except p 2 and p 3 [p~7|| . 

In the non-relativistic case the expression for a simplifies further, and it becomes possible to do one more integration. 
In this case on gets 

± j 2m(p 3 -m^±p 2 (2m-p 3 ) for annihilation, 
I Pi+p 3 ~ 2 P2 =i=2 P2 f or sca ttering. 

In the case where the light particles are in equilibrium we can further integrate over either dps (annihilation) or dp 2 
(scattering). We shall treat the case of annihilation first. For the case of a constant matrix element one gets the result 

\j(puP2,P3)t^l ] (hh ~ A/a) = G 2 ^Fi(Pi,P2), (A7) 



where 

Fi(pi, P2 ) = -/(pi)/(pa) + e-^ +E ^/ T . (A8) 
For a matrix element that is proportional to v 2 cl , one gets instead 



2 1 J2 
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\j{PuP2^)t^LiUhh A/a) = (A9) 



Using these expressions we have reduced the pair-annihilation integral to a one-dimensional integral which is quite 
straightforward to calculate numerically, 

Cann = y-^— [ P^PllFAp u p 2 ). (A10) 
(2tt) 3 2E 1 J 2E 2 1 ^ i ' I> ^ > 

Also, it is easy to see that the velocity dependent matrix element is suppressed relative to the constant one by a factor 
p 2 /to 2 , which is large in the non-relativistic case. From the above equations one can also easily calculate the velocity 
averaged annihilation cross sections for the two cases if one assumes scattering equilibrium at all times 

H»l>=G>x{=# (All) 

These can be plugged into the standard integrated Boltzmann equation 

n + 3Hn = -{(j\v\){n 2 - n 2 q ). (A12) 

Finally, the evaluation of the scattering integrals in the non-relativistic limit is also quite straightforward. Here, 
the integral must be over dp 2 to avoid integrating over the distribution function of the heavy particle. The result of 
this integration is 



fP2 max , 

H = / lJ(PuP2,P 3 t n ^Lll ] (f3h - A/a) = G 2 

J U2.mii. 



1 inf[l,a+] !ff _ t f ^ _ n 2 2 

'P2,mi„ V\PZ 

— exp 



exp 



\P\ -Pz\ 



2 

Pi + P3 " 



F 2 {pi,p 3 ), (A13) 
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where 

F 2 (pi,p 3 ) = -f(pi) + f(P3)e' E '- E3 . (A14) 
Altogether then we can write the elastic scattering term as 

* -jSFik/'^ (A15) 
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